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The Mass-Radius relation for Neutron Stars in f{R) gravity 


Salvatore Capozziello^’^’^*, Mariafelicia De Laurentis^’®’®’^^, Ruben Farinelli^^, Sergei D. Odintsov®’®’®^ 

^ Dipartimento di Fisica, Universitd di Napoli “Federico U”, 

Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, 1-80126, Napoli, Italy 
^ INFN Sezione di Napoli, Compl. Univ. di Monte S. Angelo, Edifieio G, Via Cinthia, 1-80126, Napoli, Italy 
^ Gran Sasso Science Institute (INFN), Via F. Crispi 1, 1-67100, L’ Aquila, Italy 
'^Institut fur Theoretische Physik, Goethe-Universitat, 

Max-von-Laue-Str. 1, 60438 Frankfurt, Germany 
^ Tomsk State Pedagogical University, ul. Kievskaya, 60, 634061 Tomsk, Russia 
^Lah.Theor.Cosmology,Tomsk State University of Control Systems and Radioelectronics(TUSUR), 634050 Tomsk, Russia 
^ISDC Data Center for Astrophysics, Universite de Geneve, chemin d’Eeogia 16, 1290 Versoix, Switzerland 
^ Institucid Catalana de Recerca i Estudis Avangats (ICREA), Barcelona, Spain, and 
^Institut de Ciencies de I’Espai (lEEC-CSIC), Campus UAB, 

Torre C5-Par-2a pi, E-08193 Bellaterra, Barcelona, Spain 
(Dated: December 2, 2015) 

We discuss the Mass - Radius diagram for static neutron star models obtained by the numerical 
solution of modihed Tolman-Oppenheimer-Volkoff equations in f{R) gravity where the Lagrangians 
f{R) = R + aR^{l + yR) and /(R) = R^+^ are adopted. Unlike the case of the perturbative 
approach previously reported, the solutions are constrained by the presence of an extra degree of 
freedom, coming from the trace of the field equations. In particular, the stiffness of the equation of 
state determines an upper limit on the central density pc above which the the positivity condition of 
energy-matter tensor trace T™ = p — 3p holds. In the case of quadratic f(R)-gravity, we find higher 
masses and radii at lower central densities with an inversion of the behavior around a pivoting pc 
which depends on the choice of the equation of state. When considering the cubic corrections, we 
hnd solutions converging to the required asymptotic behavior of flat metric only for 7 < 0. A similar 
analysis is performed for /(R) = R^^' considering e as the leading parameter. We work strictly in 
the Jordan frame in order to consider matter minimally coupled with respect to geometry. This fact 
allows us to avoid ambiguities that could emerge in adopting the Einstein frame. 

PACS numbers: 98.80.-k, 04.50.Kd 


I. INTRODUCTION 

The structure of a neutron star (NS) is strictly corre¬ 
lated with the equation of state (EoS), i.e. the relation 
between pressure and density in its interior j3|. Given 
an EoS, a mass-radius (M — TZ) relation and a corre¬ 
sponding maximal mass can be derived, in principle, for 
any NS. Furthermore the knowledge of these parameters 
provide significant information related to the mechanism 
responsible for NS formation and possible effects on the 
evolutionary history of their progenitors. For an intro¬ 
duction to the theory of relativistic stars, see for example 

i- 

Up to now, the physical properties of matter in strong 
gravity regimes can be uniquely studied considering the¬ 
oretical models since it is not possible to produce sim¬ 
ilar environments in laboratory. Due to this situation, 
there are more than 100 candidates for EoS, but only 
some them should be reliable once observational probes 
fix them. Consequently, the measurement of NS masses 
is thus important for our understanding on the matter 
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behavior in extreme regimes. Beside these considera¬ 
tions, it is well known that Chandrasekhar, considering 
degenerate matter, fixed a theoretical upper limit for the 
stability of a non-rotating NS at I.MMq Q. From an ob¬ 
servational point of view, the determination of mass can 
be achieved with accuracy only for NS in binary systems. 
In particular, the most accurate mass measurements have 
been derived for the binary radio pulsars where values of 
masses are around 1.35M0 [3|- However, for the X-ray 
pulsar Vela X-1, it has been measured a a mass of the 
order 1.86 ± 0.16A^ [E0 s-iid, for 4U 1822-371, a mass 
of the order 2 M 0 Q. More massive NS were discovered, 
such as the millisecond radio pulsar J0751 -I- 1807 with 
a measured mass of 2.1 ± 0.2Mq Furthermore, the 
pulsar PSR J1614 — 2230 0 has set rigid constraints on 
various EoS at strong density regimes. In summary, the 
mass of a NS cannot exceed the maximal mass limit in 
the range 3.2 T 3.6Mq according to General Relativity 
(GR) [lOL [lll |. It is important to stress that these lim¬ 
its on maximal mass exclude many EoS according to the 
observational data. Therefore, since different assump¬ 
tions provided different results, we can say that the ac¬ 
tual mass limit for NS is still a mystery. We need a more 
accurate derivation consistent with the observations. In 
other words, this state of art allows us to assume that 
the NS mass should be in the range 1.4M0 to 6Mq. This 
situation is very unsatisfactory and so there is a severe 
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need for reliable methods to obtain the NS mass limit, 
otherwise waiting for more precise observations. 

In this paper we will face the problem of NS mass 
limit adopting f{R) gravity. This is a straightforward ex¬ 
tension of GR where one relaxes the strict request that 
the gravity action is linear in the Ricci scalar R as in 
the Hilbert-Einstein case. These models can be seen as 
simple cases of a more general class of Extended Theo¬ 
ries of Gravity [l^ - [T^ . The reason why we adopt such 
an approach is that higher order curvature corrections 
can emerge in the extreme gravity regimes inside a NS 
[IS HH, US US ■ The effective pressure related to the cur¬ 
vature could naturally cure, in principle, some shortcom¬ 
ings of NS theory that are often addressed by asking for 
exotic EoS [ISl. The philosophy of the present paper is 
to construct reliable M — TZ relations solving directly the 
full modified Tolman-Oppenheimer-Volkoff (TOV) sys¬ 
tem equations [T7| . In other words, in our numerical 
integrations, we are not imposing arbitrary perturbation 
methods but solve numerically the full TOV system. In 
particular, we take into account quadratic and cubic cor¬ 
rections to the Ricci scalar and, in general, power-law 
f{R) models. The aim is to control precisely how the 
results deviate from those of GR in order to see where 
strong gravity regimes emerge and affect the M — 7?, re¬ 
lation. 

The paper is organized as follows: In Section im we 
derive the modified TOV equations for f{R) gravity. In 
particular, we consider f{R) models with quadratic and 
cubic corrections and generic power law f{R) models. 
Solutions of stellar structure equations are given in Sec¬ 
tion uni In Section lIVi we discuss the results focusing, 
in particular, on the M — TZ relation. Conclusions and 
discussion are reported in Section [V] Through the paper, 
we will indicate with TZ the radius of the object and with 
R the Ricci curvature scalar. 


(-I-, —, —, —The metric for systems with spherical sym¬ 
metry has the usual form 


ds^ = - e^^dr^ - r^{de'^ + sin^ Odcj)^), (3) 

where w and A are functions depending only on the radial 
coordinate r. Within the star, matter is described as a 
perfect fluid, whose energy-momentum tensor is T^i, = 
diag(e^’"pc^, r^p, r^psin^ 9), where p is the matter 
density and p is the pressure . The equations for the 
stellar configuration are obtained adding the condition 
of hydrostatic equilibrium which can be derived from the 
contracted Bianchi identities 


= 0, 


( 4 ) 


that give the Euler conservation equation 


dp 

dr 


-{P + P) 


dw 

dr 


( 5 ) 


From the metric m and the field equations m , it is 
possible to derive the equations for the functions A and 
w in the form 


dr 


Se^^Gnrp 


„2A 




2r{2^+rR'^) 


if 

dR 


rR" 


n/2^ 

dR3 


( 2 * + ’•R'Si 


( 6 ) 


and 


II. TOLMAN-OPPENHEIMER-VOLKOV 
EQUATIONS FOR f(R) GRAVITY 


dw Se'^^GPirr 



Let us start from the f{R) action given by 

g2A 

/-2 + (2-r2i?)^ 

-2(^ + 2rR'|^) 



2^ 


4 p 

A = ^ [f{R) + , (1) 


where g is the determinant of the metric g^i, and £matter 
is the standard perfect fluid matter Lagrangian. The 
variation of .O with respect to g^i, gives the field equa¬ 
tions dMi: 


rv. V. _ r..df{R) SttG^ 
dR 2 W/iVjy g^^vG\\ Tfiv, 

( 2 ) 

X, rj. -2 

where 1 av = j- is the energy momen- 

turn tensor of matter. Here we adopt the signature 


respectively. In both Eqs. ([51) and ([Tj), the prime denotes 
a derivative with respect to r for the Ricci scalar R{r). 

The above equations are the modified TOV equations 
that, for f(R) = R, reduce to the standard TOV equa¬ 
tions of GR (^ . [1H| . It is important to stress that, in 
/(i?) gravity, the Ricci scalar is a dynamical variable and 
then we need a further equation to solve the system of 
equations ([3), ([Hj) and ([7|). For this aim one needs to 
consider the trace of the field Eqs. m which takes the 
form 


30 


df{R) 

dR 


df{R) 

dR 


R-2f{R) 


SttG 


(P - 3p), 


( 8 ) 
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where 


is the d’Alembert operator in curved spacetime. It can 
be easily checked that for f{R) = R, Eq. ([5]) reduces to 
the trace of GR, i.e. 


R = - — {p-ip). (10) 

In order to close the system of Eqs. ([S]), ([HI), © and I©, 
one needs to provide an EoS, P{p) relating the pressure 
and the density inside the star. 

In the next subsection we will take into account two 
physically relevant f{R) Lagrangians with the aim to 
obtain the M — TZ diagram for NS. We use a fully self- 
consistent non-perturbative approach to solve the stellar- 
structure equation in their exact form. 


A. The case f{R) — R + aR^{l -|- jR) 

We first consider here a quadratic form of f{R) with cubic 
corrections, according to 


f{R) = R + aR^{l + jR), (11) 

where both a and 7 have dimensions cm“^. This class 
of models can be related to the presence of strong gravi¬ 
tational fields and emerge, for example, in cosmology, to 
achieve inflation. In particular higher-derivative curva¬ 
ture terms naturally lead to the cosmic accelerated ex¬ 
pansion of inflation ( 3 ^ . At fundamental level, their ori¬ 
gin is related to the effective actions coming from quan¬ 
tum gravity or from quantum field theory formu¬ 
lated in curved spacetime |34| . They lead to renormaliz- 
able models at the one-loop level. In the extreme grav¬ 
ity regime of NS, also if very far from the full quan¬ 
tum gravity regime, it is realistic to suppose the emer¬ 
gence of curvature corrections to improve the pressure 
effects. It is important to stress that such a model can¬ 
not be confronted with the Solar System tests of GR since 
the quadratic and cubic terms emerge in strong gravity 
regime as discussed above. They cannot be present at 
Solar System scales since the only relevant term in the 
weak field regime is the linear Ricci scalar R. The in¬ 
terest of this model, in the present context, is that the 
interior of a NS is a natural laboratory where high curva¬ 
ture regimes can emerge and lead the M — TZ relation (see 
also the discussion in [^). In some sense, the interior of 
a NS is similar to the conditions of the early universe. 

For this model, Eqs © and © take the form 


dX Ae’^^Girrp 

~ c2 [1 -h Ra{2 + 3i?7) -f ra(l -f 3Rj)R'] 

1 

Ar [I + Ra{2 + 3 R 7 ) -f ra(l -f 3 ^ 7 )^'] ^ 

{e^\Ra {R (r^ - 67 -f 2r'^R'j) - 4) - 2] 

+2[l + ZR^a -1 -f 2ra {R' (2 -f Zr-^R') + rR")] (12) 

+ARa[l + 3 r 7 (2i?' -f rR")]}, 

and 

dw 

~ [c4(l -h Ra{2 + 3 R 7 ) -f ra(l -f ?,R-^)R'] ^ ' 

1 

'^4r-[l -f Ra(2 + 3Rj) + ra(l -f 3 ^ 7 )^'] ^ 
e^^l2 + Ra(4 - R(r^ - 67 -f 2r^R'y))] 

-2[1 -f Ra(2 + 3Rj) + 4ra(l + 3Rj)R'], 

while the trace Eq. © becomes 


d^R _ 1 

d ^2 6c^ra(l -f 3 R 7 ) 
e^^r[c^i?(i?^a7 — 1) -f 8Gtt{—3P + c^p)] 
-6c^aR'[3r-fR' -f (1 -f 3R-f){2 + rw' - rA')] 


(14) 


B. The case f{R) = ^1+" 

Another interesting class is 

/(i?)=i?i+% (15) 

that are power-law models. If we assume small deviation 
with respect to GR, that is |e| <g; 1, it is possible to write 
a first-order Taylor expansion as 

R^+^ ~ R + eR\ogR + 0{e^), (16) 

which is interesting in order to define the right physical 
dimensions of the coupling constant and to control the 
magnitude of the corrections with respect to the standard 
Einstein gravity. A Lagrangian form like that in Eq. (fTBl) 
has been widely tested giving interesting results starting 
from Solar System up to cosmological scales. Applica¬ 
tions have been found in the case of the cosmological 
background of gravitational waves (^ . or in comparing 
the effects of small deviations on the apsidal motion of a 
sample of eccentric eclipsing detached binary stars (^ . 
Very strong bounds have been worked out from null and 
timelike geodesics in the cases of Solar System • Fur¬ 
thermore, black holes solutions using also Noether sym¬ 
metry approach have been found for this kind of La¬ 
grangian [ 3 ^ . [ 3 ^ . The interest of these models is related 
to the fact that the value of the parameter e can straight¬ 
forwardly relate a weak field curvature regime (e ~ 0) to 
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a regime where strong curvature effects start to become 
relevant (e ^ 0). In fact, as shown in [s^l, the Solar 
System constraints give essentially e —>■ 0. This is not 
the case for NS where high curvature regimes are well far 
from the Solar System weak field limits. In this perspec¬ 
tive, e could be different from zero in NS and then probe 
deviations with respect to GR. In this sense, NS could 
constitute a formidable test for alternative gravity. 

The explicit form of Eq. ([5]) and Q for the action (fTBl) 
takes the form 


d\ Se^^G'KrRp 

dr [2i?(l + e + elogi?) -I- reR'] 

1 

2ri?[2i?(l -I- e -I- elogi?) -I- reR'] 

{e'^^R'^lr'^Re - 2{1 + e + elogR)] 

+2[R^{1 -h e -h elogi?) - r^ei?'^ + rRe{2R' + ri?")]}, 


and 


dw _ Se'^^GPnrR 

dr c'^ [2i?(l -I- e -I- elogi?) -I- rei?'] 

1 

2r[2i?(l -I- e -I- elogi?) -I- rei?'] 
{e^^i?[r2i?e - 2{l + e + elogi?)] 
-|-2[i?(l -I- e -I- elogi?) -I- 2rei?']}, 

while the trace equation is given by 


d^R 


- + R' 

i?[c4i?(l - e) 



-k 8G7r(3P - 


3c^e 


(19) 

c^p) + c^i?elogi?] 


In the framework of these models, let us now discuss the 
stellar structure for this models with the aim to achieve 
the M — TZ relation for some physically relevant EoS. 


III. THE STELLAR STRUCTURE EQUATIONS 


Let us consider dimensionless variables for solving the 
system of Eqs. (O-dH]). We set the dehnitions x = r/rg, 
R = R/rg^, p = P/Pq, p = p/po, where = GMq/c^, 
Pq = Mqc^/ r^ and po = where = 1.48 x 10® 

cm. By the substitution R' = q, the trace equation for 
i? is lowered by an order of derivation, and the full set 
of equations is reduced to a system of first-order ODEs, 
which can be expressed as 


X' = Fi{X,w,'w',R,q,q',p,x), (20) 

w' = F 2 (A, X',w',R,q,q',p,x), 

R' = F3{q,x), 
q' = F 4 {X,X',w,w',R,q,p,x), 
p' = F5 {w',p,x). 


In such a way, we can deal with the stellar structure 
equations under the standard of dynamical systems (see 
also [2l|). 

The requirement of asymptotic flatness for the met¬ 
ric implies that A —>■ 0, u> —>■ 0, i? —^ 0 for r —>■ oo. The 
boundary value problem (BVP) can be reduced to an ini¬ 
tial value problem (IVP) by imposing initial conditions 
at the star center a; = 0 such that the function asymptot¬ 
ically have the needed behaviour. For the pressure, the 
condition is simply p(0) = Pc, where Pc is determined by 
the chosen EoS once the central density Pc is hxed. For 
the remaining functions, the problems deserve some con¬ 
siderations: first of all, we notice that the natural initial 
conditions for A and the Ricci scalar are A(0) = 0 and 
i?'(0) = 0. On the other hand, it is worth pointing out 
that the potential function w appears in Eqs ©-(HI) only 
through its first and second derivatives w' and w". This 
means that its initial value ui(0) can be defined up to an 
arbitrary constant. Once the full solution is obtained, the 
value of w at the star center can be shifted a posteriori 
such that w —>■ 0 asymptotically. 

The most critical point for the numerical solutions con¬ 
cerns the value of the Ricci scalar at the center. In 
the classical version of the shooting method, the right 
initial value of Rq is given by the root of the function 
F (,Ro) — 2 /(?? 0 ;^max) R (^max); where y(?? 07 ^max) is 
the value of i? at the boundary of the domain obtained 
with the initial condition i? 0 j while i?®(a:niax) is the im¬ 
posed boundary condition (BC). Once provided an initial 
guess for Rq, algorithms for root-searching, such as the 
Newton method for quadratic convergence, give the ini¬ 
tial values, if they are not too far from the solutions. 

The implicit assumptions of the shooting method how¬ 
ever is that the function y{Ro, a^max), and in turn F{Rq), 
is defined for any Rq ; moreover the classical computation 
requires to determine numerically the derivative dy/dRQ 
at each iteration. 

In our case, this approach has not been possible for 
two reasons: first of all, for some initial guess values 
of Rq, smaller than the right one, the function R at 
some distance x > Xcrit becomes negative, leading to 
unphysical solutions for the system of Eqs ([5|)-([5]). Ad¬ 
ditionally, the numerical computation of the derivative 
dy/dRQ at a;niax becomes time-consuming and strongly 
ill-conditioned, depending critically on the choice of the 
step-size for Rq. The latter problem originates because 
of the diverging behavior of the negative i?-function at 
values progressively higher than icrit- 
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We thus implemented a different method as follows: we 
first define for Rq an initial guess interval ~ i?o/5 
and ft! 5i?o, where Rq = 87r(pc — 3pc) is the GR 

value. Then we use a bisection method over the given 
interval, which allows to get progressively closer to the 
right value of Rq through subsequent interval subdivi¬ 
sions. The iteration stops when a further interval subdi¬ 
vision does not produce appreciable change at i?(a:max)- 
The needed accuracy in the determination of Rq depends 
on the form of /(i?) and its variable parameters. In the 
case of /(i?) = R + aR^{l + yi?), we found that, for all 
considered values of a, the fine-tuning needs to be very 
tight, independently of the value of 7 . 

However, the bisection procedure over the initial inter¬ 
val of i?o-values allows to achieve a precision as high as 
desired, which is in fact limited only by the numerical ac¬ 
curacy of the integration routines. In general, a few tens 
of iterations are required leading to the determination of 
Rq with a precision of order ARq/Rq < 10“^. For the 
Lagrangian f{R) = R^~^’^, on the other hand, as e gets 
progressively closer to zero, the results become much less 
sensitive to Rq; namely, Rq = with k of the or¬ 

der 2, gives the same stellar masses and radii. In this 
case we considered the value Rq = l/2(i?™™ -|- 
The adimensional NS radius is defined by the condition 
pixns) = 0 ; for X > Xns the integration is continued in 
vacuum for the metric potentials A and w, and for the 
Ricci scalar R. 

To check the accuracy of the results, we tested both a 
Gear M = 2 method and 4*''-order Runge-Kutta method 
with adaptive step-size and accuracy control of the func¬ 
tions during the integration of the set of above differential 
equations (see [13 for details); we found no significative 
differences between the two methods. Finally, the gravi¬ 
tational mass is obtained a-posteriori by using a Gauss- 
Legendre quadrature rule for the integral 

P^ns 

Mgrav = PoTg / A7Tx‘^p{x)dx, (21) 

where Xns is the adimensional NS radius, while p{x) is 
obtained by computing the inverse of the function p(p) 
of the EoS. 


IV. RESULTS 

As previously mentioned, the stellar structure equa¬ 
tions require to specify the EoS relating the matter pres¬ 
sure and density. From a mathematical point of view, the 
EoS is the algebraic relation needed to close the system. 
There is however an important point to keep in mind 
that has been surprisingly little discussed in literature, 
namely the maximal density pc at the center of the star 
which can be used in numerical simulations. 

Depending on the stiffness of the EoS, the maximum 
implementable value of pc is given by the need to preserve 
the physical condition Tr{T) = p^c^ — Sfy, > 0 for the 
trace equation. In the standard TOV system or in the 


perturbative approach, Tr{T) does not appear explicitly 
in the equations of the stellar structure. However in the 
exact treatment, Tr(T) represents the source term in the 
right hand side of the trace equation differential operator 
(see Eq. ®). Values of pc violating its positivity lead 
to un-physical results. As a consequence, the top-left ex¬ 
tension of the branches oi M — TZ diagrams is limited by 
the maximum achievable value of pc preserving the pos¬ 
itivity of the energy-matter tensor for a given EoS. This 
limits of course the possibility to investigate the stabil¬ 
ity of compact stellar structures up to arbitrary central 
density values. In this paper, we used the analytical rep¬ 
resentation of EoS with different stiffness as reported in 
[H (labeled as BSkI9,BSk20 and BSk2I) and the Sly 
EoS reported in (d^ . The maximum allowed central den¬ 
sities in units of 10^® gr/cm^ are pc = 2,1.25,1.35,1.7 for 
EoS BSkl9,BSk20, BSk21 and Sly, respectively. 

A. Results for f{R) = R-\- aR^ 

We first considered quadratic f{R) gravity as reported 
in Eq. m with 7 = 0 and progressively increasing val¬ 
ues of a. Note it has to be a < 0 because of the adopted 
sign convection of the metric signature (see Eq. ([3])). 
Adopting negative a, we are avoiding ghost modes and 
instability behaviors An example of the radial be¬ 
havior of the metric potentials A and w, the Ricci scalar 
R and the pressure P is shown in Fig. [TJ 



R(km) 


Figure 1: Radial profile for the metric potentials A and |w|, 
the dimensionless Ricci scalar R/vg and the dimensionless 
pressure P/Pq, obtained from the solution of the system of 
equations ([5|, (I12II . (I13II . and (I14|l for f{R) — R + aR^ with 
|q| = 1. The central densities reported are pc = 1 x 10^® 
gr/cm® {thick line) and pc = 1.5 x 10^® gr/cm® {dashed line), 
for the Eos BSkl9. 

The M — TZ diagrams obtained for the four different 
EoS here considered are instead shown in Eig. [21 to¬ 
gether with the solution of the classical TOV {a = 0). 
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Independently of the EoS, as |a| increases, the following 
behavior can be observed: the mass at lower central den¬ 
sity are higher, then the branch crosses the classical GR 
solution and for pc > ^f(R) < -^GR- The value 

of depends on the EoS, and it is approximatively 
0.6 X 10^® < < 1 X 10^^ (Fig. O. The results re¬ 

ported in Fig. [5] and Fig. [3] are obtained for values up 
to |a| = 20. Further higher values produce additional 
counter clockwise rotation of the M — TZ traces, until 
getting a saturation for |a| > 100. There is an anti¬ 
correlation between a and the required value of the Ricci 
scalar i?o at the NS centre simultaneously satisfying the 
asymptotic condition i? —> 0. This can be qualitatively 
understood also looking at the trace of the field equations 
in quadratic /(i?)-gravity, which can be written as 

UR + ni^{R + xT) = ( 22 ) 

, 2 1 ■ ^ . SttG , . 

where m =-is an euective mass, v = —and, m 

6 a . . ’ ^ c4 

vacuum, the i?-function exponentially decays approxima¬ 
tively as 6 “'"''. The higher the value of |a|, the higher 
is the distance at which the Ricci scalar asymptotically 
goes to zero. A function obeying the condition i?'(0) = 0 
and the asymptotic zero-value condition at large radii 
needs to be smaller at r = 0 when the exponential-decay 
constant m decreases. 

B. Results for f{R) = R + aR?{l + yi?) 

Let us improve now the above considerations by taking 
into account f{R) models with non-zero cubic correction 
term with 7 7 ^ 0 in the Lagrangian (fTTl) . It is also ev¬ 
ident that the higher the value of |a|, the higher must 
be the value of jyj in order to see appreciable deviations 
from quadratic /(i?)-gravity. We discuss here the case 
|a| = 0.5 as an example, showing the results in Fig. 0] 
for the case of a SLy EoS. For positive and progressively 
increasing values of 7 , the behavior of the traces in the 
M — TZ diagram is similar to that of the quadratic f(R)- 
form at increasing values of |a| - higher masses and radii 
below a certain central density pc and inversion of the 
trend above it. This results can be understood by keep¬ 
ing in mind the definition of f{R) in Eq. (fTTll . where 
the cubic term a^R^ is negative for 7 > 0 (being a < 0 
because of the adopted signature). 

Thus, positive 7 -values provide a further decrease of 
the value of the Lagrangian density f{R) in the same 
fashion as the quadratic coefficient |a| does. This nega¬ 
tive contribution unavoidably reflects in the M — 7?. re¬ 
lation both for the quadratic and quadratic plus cubic 
corrections in R. 

At the opposite, negative 7 -values produce a clockwise 
rotation of the M — TZ trace. This can be seen indeed in 
Fig. [Ufor the case 7 = —10. We found however that if 
from one side 7 can take arbitrary high negative values, 
this is not true in the opposite case. For instance, for 
the case |a| = 0.5 and 7 = 20, here discussed, we were 


not able to find converging solutions matching the re¬ 
quired asymptotic behavior. As outlined above, the crit¬ 
ical maximum value of positive 7-values depends however 
on the simultaneous choice of a. As an example, the com¬ 
bination |a| = 57 = 20 gives rise to a converging solution 
which is however almost indistinguishable from the case 
with 7 = 0. 

In general, we can say that for each value of |a|, there 
is a critical positive value above which solutions 

are not allowed. More specifically, for 0 > 7 > 7'^''** 
quadratic and cubic forms of /(i?) provide very similar 
results, while for 7 < 0 it is possible to go beyond val¬ 
ues where modifications of the M — TZ relation are more 
evident. In the perturbative approach both positive and 
negative values of 7 were considered 0- 

C. Results for f{R) = R + eR\ogR 

Finally, we report the results for the model given in 
Eq. (flSll . Albeit the field equations given in Sec. IllBl 
are presented for clarity using the Taylor expansion for 
f{R) of Eq. (II6D . it is worth noticing that we actu¬ 
ally solved the stellar structure equations using the exact 
form f{R) = R^'^^ . Similarly to Fig. [T 1 we show, in Fig. 
[SJ an example of the radial behavior of the functions A, 
w, R and P. The reason why we considered the case 
e <C 1 can be better understood looking at the M — TZ 
diagram presented in Fig. |6l For |e| < 0 .01 a significant 
deviation is observed with respect to the classical TOV; 
at first sight, it is noticeably that the traces in the dia¬ 
gram present a self-similar behavior for increasing values 
of |e|. The results of physical interest (higher masses and 
radii) are obtained for e < 0, while, in the case of pos¬ 
itive values of e, the traces are blended with respect to 
the classical TOV solution. Note that in this case there 
is no pivoting around some critical central density, unlike 
what found in the case of the polynomial form of f{R). 
Actually the resizing of the action f{R) = reflects 

in a resizing of the traces in the M — TZ diagrams. In 
the case of Eos BSK 20 and BSK 21 , NS masses, larger 
than3M0, can be easily achieved; these values must also 
be considered as lower limits on the NS mass, which can 
be even higher for fast rotating objects. 


V. DISCUSSION AND CONCLUSIONS 

Neutron Stars have a main role in relativistic astro¬ 
physics for several reasons. They are the most stable 
compact objects of the universe (a part the black holes) 
where matter reaches extremely high field regimes. Be¬ 
sides, they are very important in order to understand 
the final stages of stellar evolution. It is important to 
point out that the extreme field regimes in NS cannot be 
achieved in any ground-based lab so they have acquired 
a crucial role also in nuclear and particle physics. Due to 
this feature, understanding the EoS working into a NS 
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Figure 2: M — TZ diagrams for different EoS (labeled in the top of each panel) and quadratic form f{R) = R + aR?. For each 
EoS the maocimal central density is determined by the condition pc — 3p > 0 (see text). The classical TOV solution is also 
reported. The values of a are in modulus. 


has a twofold meaning: from one side, it can give infor¬ 
mation on the state of matter in this compact objects, on 
the other side, it is relevant at astrophysical level to un¬ 
derstand the global behavior of stars in the final stages of 
their life. Besides, these standard roles, NS could be ex¬ 
tremely important to test alternative theories of gravity, 
due to the huge gravitational field acting on them. 

The aim of this paper is to show that the M — 77. re¬ 
lation of NS can be consistently achieved by extended 
theories of gravity as f{R) gravity. In particular, we pro¬ 
posed some physically relevant f{R) models and modi¬ 
fied the TOV equations accordingly. In the stellar struc¬ 
ture equations, new terms related to curvature correc¬ 
tions come out and lead the evolution of the M — 77 rela¬ 
tion. Physically, these terms assume the role of a sort of 
curvature pressure capable of leading the mass and the 
radius of the star Specifically, we reduced the stellar 
structure equations to a dynamical system and integrated 
it numerically putting in evidence the role of curvature 
corrections into the integration. The resulting M — 77 
diagrams strictly depends on the value of the curvature 
corrections, the sign of the correction parameters, and 
the chosen EoS. We dealt with these terms as corrections 
to the GR in order to control deviations with respect to 
the standard Einstein theory. 

This point deserves a further discussion according to 


the numerical results presented in the above figures. Let 
us consider first Figl^l where the M — 77 relation is re¬ 
ported for several values of the parameter a and different 
EoS. The GR value is for a = 0 while a ^ 0 represents 
corrections with respect to GR. Glearly the increasing of 
a in modulus gives rise to a sort of stretching-blending 
rotation of the M—77 relation slope around a fixed M—77 
point that sits in the intervals (1 t1 .2)M0 and (11 t 12.5) 
km. This means that the magnitude of the gravitational 
corrections alters the structure of the NS thanks to the 
further curvature pressure-density terms present in the 
TOV equations (see also 0 ). In some sense, given a 
EoS, curvature corrections result relevant in shaping the 
stellar structure and determining the M — 77 relation. 
Specifically, the stretching-blending rotation depends on 
the effective pressure related to the curvature corrections 
in Eq. and, in particular in Eq. (fl^ . Increasing 
the parameter a means that the original M — TZ relation 
of GR is affected by a further pressure term that mod¬ 
ifies the effective mass of the star and consequently is 
radius. However, also the effective density results mod¬ 
ified by the curvature corrections. Then the net effect 
is the stretching-blending rotation around M and 77 val¬ 
ues where curvature corrections are not so relevant (see 
Figl2). It is important to stress again that the stretching¬ 
bending effect is due to the curvature corrections and not 










BSk19 BSk20 



Figure 3: Mass versus central density for different EoS (labeled in the top of each panel) and quadratic form f(R)=R+QR^. For 
each EoS the maximal central density is determined by the condition pc — 3p > 0. As in Fig. [2j the modulus of a is reported 
(see text). 


SLy 



Figure A: M — TZ diagram for f{R) = R + aR^{l + -yR) with 
|q| = 0.5 and different values of 7 , and Sly EoS. 


to the change of standard matter EoS. 

Similar considerations hold also for the M — pc dia¬ 
grams (Figl3]) where the total mass of the NS is given as 
a function of the central density pc- Also in this case we 
have a stretching-blending rotation around given values 
of mass and density (GR values) depending on the abso¬ 
lute value of a. Inserting also the cubic correction, led 



0 2 4 6 8 10 12 14 16 18 20 


R (km) 

Figure 5: Same as in Fig. [T]but for f{R) = R + eRXogR with 
£ = -0.05. 


by the parameter 7 , the trend is similar (see Figd]). 

The case of gravity is totally different. In this 

case, as we can see from Fig El different values of the 
parameter e make to scale the M — TZ relation so that 
larger stable structures, in terms of mass and radius. 
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Figure 6: M — TZ diagrams for f{R) given in equation (I15|l with four EoS and different values of e. The classical TOV 
corresponding to e = 0 is also reported. 


are achieved also if slight variation with respect to the 
value e = 0 are considered. This fact could be extremely 
relevant in order to address extreme massive NS as re¬ 
vealed by some observations Specifically, while for 
the quadratic and cubic models, the M — TZ relation is 
only modihed (stretching-blending rotation) just chang¬ 
ing the stability region of NS with respect to GR, here 
the curvature pressure and density give rise to a scaling 
law. This is quite obvious due to the model, but the 
physical consequences are relevant since extreme massive 
and large objects can be achieved, as one can see from 
Fig. [HI In such a case, the breaking of GR behavior, re¬ 
lated to e ^ 0, gives rise to a new stability branch for NS 
allowing extreme objects. We stress again that this phe¬ 
nomenon is strictly related only to the effective curvature 
quantities and not to the change of standard matter EoS. 

It is important to stress that we never used exotic mat¬ 
ter but only realistic EoS. This point is crucial since we 
adopted always the Jordan frame so that the gravity sec¬ 
tor results corrected while the matter sector is unaffected. 
In such a way, the geodesic structure is not altered and 
standard EoS can be assumed. On the contrary, if we 
were adopting the Einstein, frame, we would have the 
standard gravity sector but a non-minimally matter sec¬ 
tor. In such a case, geodesic structure results altered and 
it could be dangerous to adopt standard EoS. In other 
words, as it is shown in detail in (d^ , the conformal trans¬ 


formations have the effect of shifting the non-minimal 
coupling (in our case f'{R)~^) from the gravitational 
to the matter sector and then the meaning of quantities 
like gravitational potentials, pressure, matter density and 
mass result more complicated and have to be accurately 
discussed. In summary, we develop our calculations in 
the Jordan frame considering it as the "physical" frame 
and avoiding the ambiguities that could emerge in the 
Einstein frame (see, for ex amp le, (d^ l^l. The same 
approach is adopted also in . 


As an example of these considerations, we can see that, 
qualitatively, the behavior of the traces in the M — TZ dia¬ 
grams with respect to GR is the opposite with respect to 
that reported by [d^ for the same form f{R) = R + aR^ 
(see FiglHand Figs.l and 2 in (d§|). The only possibil¬ 
ity to explain this crucial differences relies on the fact 
that computations in were performed in the Einstein 
frame, while in our paper we work in the Jordan frame. 
The only way to compare exactly the results is to com¬ 
pare the behavior of the M — TZ diagram under conformal 
transformations assuming the same EoS. This will be the 
topic of a next project. Finally, the results of this work 
will be extended to magnetic as well as rotating neutron 
stars. 






adopted the analytical form 


10 


Appendix A; Analytical representations of 
Equations of State 


Here we report the functional form of EoS that we used 
along the paper to solve numerically the stellar structure 
equations. The pressure can be parameterized as a func¬ 
tion of density. Let us denote with ^ = log(p/g cm“^) 
the dimensionless density and with ^ = log{P/dyn cm“^) 
the dimensionless pressure. We used the SLy equation as 
reported in for non-rotating NS configurations. It is 

oi -ha2^-f 03^3 

C = -j-Jo a5 U-«6 

+ {ar + asO Magiaio - 5)) 

-|-(aii -I- 0125) /o(ai3(oi4 — 5)) 

-|-(ai5-I-0165)/o(ai7(oi8 — 5)) j (Al) 

where the parameters oi for SLy EoSs are given in Ta¬ 
ble U For BSkl9, BSk20 and BSk21, following [Hi, we 


01 - 1-025 + 035 ^ 1 

C = -r--?- {exp 05(5 - Oe) + 1} 

1 + 045 

+ (07 + 085 ) (exp [09(06 - 5 )] + 1 }”^ 

+ (oio + oii5) (exp [012(013 - 5)] + 1}“^ 

+ (oi4 + 0155) {exp [oi6(oi7 - 5)] + 1}“^ 

^ 1 + [019 (5 - 02o)]2 ^ 1 + [022 (5 - 023 )]^^ ^ 


The values of parameters Oi are given in Table HIl 


Table I: Values of ai(SLy) parameters for Eq. (lAll) . 


i 

Hi (SLy) 

i 

Oi(SLy) 

1 

6.22 

10 

11.4950 

2 

6.121 

11 

-22.775 

3 

0.005925 

12 

1.5707 

4 

0.16326 

13 

4.3 

5 

6.48 

14 

14.08 

6 

11.4971 

15 

27.80 

7 

19.105 

16 

-1.653 

8 

0.8938 

17 

1.50 

9 

6.54 

18 

14.67 
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